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Abstract. - We present a numerical study of magnetic phases of the 2D electron gas near 
freezing. The calculations are performed by diffusion Monte Carlo in the fixed node approxima- 
tion. At variance with the 3D case we find no evidence for the stability of a partially polarized 
phase. With plane wave nodes in the trial function, the polarization transition takes place at 
r 3 = 20, whereas the best available estimates locate Wigner crystallization around r s — 35. 
Using an improved nodal structure, featuring optimized backflow correlations, we confirm the 
existence of a stability range for the polarized phase, although somewhat shrunk, at densities 
achievable nowadays in 2 dimensional hole gases in semiconductor heterostructures . The spin 
susceptibility of the unpolarized phase at the magnetic transition is approximately 30 times 
the Pauli susceptibility. 



Quantum simulations show that the three-dimensional electron gas undergoes a continuous 
spin-polarization transition, upon decreasing the density into the strongly correlated regime, 
before forming a Wigner crystal jjl p|. Despite the theoretical interest of a simple model 
exhibiting quantum phase transitions |0] , the difficulty of realizing a low density electron gas in 
real materials makes the contact with experiment a rather indirect one Q. However, electrons 
can be confined into effectively two-dimensional systems, for instance in Si MOSFET's and 
III-V semiconductor heterostructures Q , over a density range extending down to the freezing 
transition ||. 

Far from being a mere playground for testing many-body theories and numerical simu- 
lations, strongly correlated two-dimensional electronic systems offer an extremely rich and 
interesting phenomenology, such as the fractional quantum Hall effect [[?J and a previously 
unexpected metal-insulator transition not to mention their relevance to superconductor 
cuprates S. Therefore, we think that it is useful to extend the work of ref. and study the 
polarization transition of the two-dimensional electron gas (2DEG). In fact, spin fluctuations 
appears to play an important role in the 2DEG at the metal insulator transition in presence 



of disorder |10|) 
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Fig. 1 - Optimized pseudopotentials it M „(r) (solid lines, left scale) and VMC pair distribution functions 
9nv( r ) (dashed line, right scale) for the unpolarized 2DEG at r s = 20. 



Previous quantum Monte Carlo (QMC) simulations JTl]— [l5|] of the 2DEG provide little in- 
formation on the magnetic properties of the system, with only one spin susceptibility calcula- 
tion p5[ at non-zero wave vectors and relatively high densities (r s < 10, where r s as = l/^/np 
with as the Bohr radius and p the density), and one attempt at exploring the partially polar- 
ized case near freezing. We emphasize that we have only mentioned, here, results obtained 
from QMC, a computationally demanding numerical method which projects the ground state 
from a trial function by a stochastic technique jl6j . Both in this work and in refs. [p|,pl]^15[| 
the method is implemented in the fixed-node (FN) approximation , giving the lowest up- 
per bound to the ground state energy consistent with the nodal structure of Wy. Though not 
exact unless \I/t has the same nodes as the unknown ground state, a condition not generally 
met, FN-QMC provides fairly accurate predictions for ground state properties of interacting 
fermion systems. Other approximate theories may fulfill such an accuracy requirement at 
weak coupling |Q, but are not amenable to systematic improvement in the strongly corre- 
lated regime. When applied beyond their range of validity, they predict a variety of magnetic 
instabilities at unrealistically high densities [BOj . 

We simulate 2D systems of N = Nf + 7Vj, electrons, interacting with a pair potential 
v(r) = 2/(r s r), in a square box of side L = VjVtt with periodic boundary conditions (energies 
are in Ry and distances in r s as units). A rigid uniform background of positive charge ensures 
neutrality. The polarization is £ = (JVj — iVjJ/iV. Following ref. S, we use a Jastrow- 
Slater form for the trial function, — exp[— J2i<j u ^( r ij)]^T^l- Here the indices i,j label 
electrons, /i and v label the spin of particles i and j; u(r)^ v are pair pseudopotentials for like 
(/i = v) or unlike-spin (/i ^ v) particles, optimized for each system considered, and -D-j-q) is 
a Slater determinant of plane waves for up(down)-spin electrons. As it has been previously 
done in three dimension we are only considering homogeneous phases fLgf| . 

Figure |l| shows the optimized pseudopotentials for an unpolarized system with TV = 114 
at r s = 20. In the same figure, the spin-resolved pair distribution functions g-\-\{r) and g-\i{r), 
computed with variational Monte Carlo (i.e. by taking the expectation value over \^t)) is 
shown by the dashed lines. Despite the fact that the like-spin pseudopotential is less repulsive, 



D. Varsano, S. Moroni and G. Senatore : Spin-Polarization transition in etc. 3 




-92.9 



-63.8 - 



g -63.9 I-' 



-64.1 



^ I I I 1 I , I 1 I , I 1 I , I , I , , 










- 1 < 

-* 

_s 


9 9 

9 


* * 


8" 

s - 


8 






s — 


1 { 


data 








1 , , 


fit 
extr 

1 , , 1 


1 , , 1 


» 





0.2 0.4 


0.6 0.8 


1 



Fig. 2 - Plane wave fixed node energy versus polarization at r s — 20 (left panel) and at r s = 30 (right 
panel). 



the pair distribution functions do not show the local ferromagnetic order predicted S in 3D 
across the magnetic transition. Note that in 2D, at the accuracy level of FN-QMC with 
plane wave determinants, the energies of the unpolarized and the fully polarized phases cross 
precisely at r s = 20, according to both ref. [O and the present work (see below). 

Our estimate of energy vs. polarization in the thermodynamic limit, parametrized in the 
form E(C) = E{C, = 0) + /?C 2 + 7C 4 , is shown by the solid line for r s = 20, 30, in fig. |. The 
statistical uncertainty is less than 10 -5 . The points with errorbars are FN-QMC results with 
plane waves, and the open circles are the result of the fit for the number extrapolation (see 
below) . The arrows indicate the results of ref. [[u} at £ = and £ = 1. 

The polarized and unpolarized phases have nearly the same energy at r s = 20, whereas 
at r s = 30 the polarized fluid is stable. We find no evidence for the stability of a partially 
polarized phase at these densities, and we believe that this is the case for other densities as well 
(a reentrant partially polarized phase could form for r s < 20; alternatively, 7 should change 
sign twice between r s — 20 and r s = 30: both cases seem unlikely). This conclusion 
agrees with a conjecture of ref. [0, based on the simulation of a single partially polarized 
system at r s =40; however, the more systematic analysis presented here seems worth the 
effort, also in view of some discrepancies existing between the results of ref. Jll]] and later 
simulations fil||l4f|. 

The dispersion of the FN E(() is very weak, which reflects in a large spin susceptibility \s 
of the unpolarized 2DEG at r s — 20, where the first-order polarization transition takes place. 
The best-fit value of (3 is 8.2±1.2xl0 -5 , corresponding to Xs/xp — 30, where xp 1S the Pauli 
susceptibility. This is an extremely large value of Xs- We note that a phenomenological model 
for Xs 01) based on including known sum rules and fitting QMC data from ref. fu]| , gives 
Xs/xp = 1-6 at r s = 20. 

We now discuss the differences between our results and the ones of refs. 13, 14| (all three 
calculations are nominally equivalent). We perform QMC simulations in the Diffusion Monte 
Carlo (DMC) version p3fl . A population of N\y walkers (copies of the system) perform a 
random walk in configuration space. The time steps of the random walk stochastically follow 
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Table I - Energy-fitting parameters from DMC. See text for the meaning of symbols 





E(C = 0) 


(3 x 10 5 


7 x 10 s 


A 


e 


b 


a 


x 2 


r s = 20 


-0.092456 


8.22 


-7.34 


0.00129 


2.0 


0.00712 


-0.0030 


1.7 


r s = 30 


-0.063783 


-0.64 


-2.72 


0.00045 


1.8 


0.00397 


-0.0007 


1.9 



a short time approximation to the imaginary time evolution of the quantum system. The 
results need to be extrapolated to zero time step r and, in the implementation used, to 
infinite Nw, as well as to the standard thermodynamic limit N — > oo. In brief, (i) we perform 
the number extrapolation directly on several DMC energies instead of sticking the VMC size 
corrections on the DMC energy computed for a single system size E^H, and (ii) we believe 
that we better control the finite r and finite Nw biases. 

We fit the energies computed at various polarizations £ and particle numbers N using 

E N (0 = E(C = 0) + /?C 2 + 7C 4 + A(l + eC 2 )(4 0) (0 - T<°>(0) - (& + a( 2 )/N, (1) 

where T^ 0) and are the kinetic energies of a non-interacting 2DEG at r s — 1 for N 
particles and in the thermodynamic limit, respectively. The fitting parameters from FN 
energies with plane wave nodes are listed in table |. Obviously, the fit to the VMC energies 
yields different values for E(( = 0), /3 and 7. Less obvious is the difference of the other 
fitting parameters (i.e. is the number extrapolation the same for DMC as for VMC?). In 
table H we list the fitting parameters obtained from VMC (103 systems at each density using 
parameter-free RPA pseudopotentials (24j] to avoid spurious effects from different levels of 
optimization at different particle numbers). This VMC-DMC comparison may be in part 
biased by the fact that we used Ewald sums |2f| for the Coulomb potential in VMC, and 
a truncated potential p(| in DMC. However we note that b, the fitting parameter directly 
related to the potential energy, is nearly the same in VMC and DMC, a significant difference 
being instead present in the kinetic energy terms. 

We turn now to the time step error and the population control bias. The results for 
r s = 20 shown in fig. are all extrapolated at infinite number of walkers and zero time step. 
For r s = 30 the extrapolation has been done in 10 cases; then a time step f and a number of 
walkers Nw have been chosen such that the residual bias is less than 5 x 10 -6 in the selected 
cases. Some of the results shown in fig. || for r s — 30 are obtained using f and Nw- 

For 58 particles at r s = 20, both ref. @ and [u| report E = -0.09248(1). However we 
find that the extrapolation at zero time step and infinite number of walkers is —0.092557(76). 
In fig. |3] the results for various numbers of walkers and time steps are shown by star symbols. 
For each r the extrapolation (linear in 1/Nyy) to infinite number of walkers is shown by the 
diamond. The extrapolation (linear in r) of the diamonds to zero time step is shown by 
the open circle. The shaded region encloses the energy values between —0.09248 ± 1, i.e. it 
represents the result of refs. |l3|,^4|. 

The filled circle indicates a result, extrapolated to r = and Nw — 00, obtained with a 

Table II - Energy-fitting parameters from VMC. See text for the meaning of symbols 
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Fig. 3 - Effects of finite time step and number of walkers on the fixed node energy for 58 particles at 
r s — 20, in the normal fluid. See text for the meaning of symbols 



different DMC code, in which walkers carry a weight for a while and branch every so often 




instead of branching at each time step. The cross is the result, extrapolated to r = 0, of yet 
another calculation with reptation QMC p8| , a single-thread algorithm (hence, one having 
no population control bias). The finite-r and finite- Nyy biases are all different in these 
calculations, but the extrapolated result is the same. This cross-check performed using three 
different algorithms gives us confidence on our extrapolated results. 

The above analysis demonstrates the importance of an extremely accurate control of the 
convergence parameters of the numerical calculations, due to the relevance of small energy 
differences on physical properties related to spin polarization in this strongly coupled system. 

We finally consider the FN approximation. In 3D, the error of the FN energies with 
plane wave nodes is presumably negligible in characterizing the continuous polarization 
transition. On the energy scale of the weak dependence of E on £ shown in figs. |[ however, 
the FN bias might be relevant. We have computed FN energies with better nodes provided by 
optimized backflow correlations [S, at r s = 20 and 30 for zero and full polarization. Where 
available, comparison between the FN energy with backflow nodes and the exact result 
supports the accuracy of the former for the 2DEG. Including backflow correlations, the FN 
energy decreases more for £ = than for £ = 1: the relative gain of the unpolarized phase 
is found to be 7.2(7) x 10~ 5 at r s = 20 and 2.0(9) x 10~ 5 at r s = 30. This makes the 
unpolarized phase stable at r s — 20, but at r s — 30 the energy of the polarized fluid remains 
lower (compare with fig. || and table |), even though by an amount close to the statistical 
uncertainty 

At r s — 30 the system is still fluid. Based on plane wave FN energies, crystallization 
occurs at a lower density, r s = 37 ± 5 and 34 ± 4 according to refs. Jll| and Q , respectively. 
In fact, backflow correlations will further stabilize the fluid phases, likely yielding a slightly 
lower crystallization density (for the polarized liquid at r s = 30 we find a tiny difference of 
1.6(4) x 10~ 5 between the plane wave and the backflow FN energy). 

Summarizing, we have studied the polarization transition in the 2DEG using the FN-DMC 
method. With plane wave nodes, the transition is first order and takes place at r s = 20, where 
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the spin susceptibility reaches 30xp. Upon improvement of the nodal structure, obtained 
including backflow correlations, the transition density approaches r s — 30, still leaving a 
stability window for the fully polarized phase between the paramagnetic fluid and the Wigner 
crystal. 

* * * 
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